In this tutorial we compare the formulation of a Navier-Stokes problem using standard assembly with mixed function spaces and block assembly.
import typing
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import numpy.typing
import petsc4py.PETSc
import ufl
import viskex
import multiphenicsx.fem
import multiphenicsx.fem.petsc
INFO:root:running build_ext INFO:root:building 'multiphenicsx_85a06d4525c3410c54a93869d159d912' extension INFO:root:creating /tmp/tmpsq_50hk9/github INFO:root:creating /tmp/tmpsq_50hk9/github/home INFO:root:creating /tmp/tmpsq_50hk9/github/home/.cache INFO:root:creating /tmp/tmpsq_50hk9/github/home/.cache/fenics INFO:root:creating /tmp/tmpsq_50hk9/usr INFO:root:creating /tmp/tmpsq_50hk9/usr/local INFO:root:creating /tmp/tmpsq_50hk9/usr/local/lib INFO:root:creating /tmp/tmpsq_50hk9/usr/local/lib/python3.10 INFO:root:creating /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages INFO:root:creating /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx INFO:root:creating /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp INFO:root:creating /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx INFO:root:creating /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/fem INFO:root:creating /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/la INFO:root:creating /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/wrappers INFO:root:x86_64-linux-gnu-gcc -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -fPIC -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp -I/usr/local/lib/python3.10/dist-packages/mpi4py/include -I/usr/local/lib/python3.10/dist-packages/petsc4py/include -I/usr/local/dolfinx-real/include -I/usr/local/slepc/include -I/usr/local/slepc/linux-gnu-real64-32/include -I/usr/local/petsc/include -I/usr/local/petsc/linux-gnu-real64-32/include -I/usr/local/include -I/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration -I/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/wrappers -I/github/home/.cache/fenics -I/usr/include/python3.10 -c /github/home/.cache/fenics/.rendered.multiphenicsx_85a06d4525c3410c54a93869d159d912.cpp -o /tmp/tmpsq_50hk9/github/home/.cache/fenics/.rendered.multiphenicsx_85a06d4525c3410c54a93869d159d912.o -std=c++11 -fvisibility=hidden -std=c++20 -DDOLFINX_VERSION=0.7.0.0 -DDOLFINX_VERSION=0.7.0.0 -DHAS_ADIOS2 -DHAS_SLEPC -DHAS_PTSCOTCH INFO:root:x86_64-linux-gnu-gcc -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -fPIC -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp -I/usr/local/lib/python3.10/dist-packages/mpi4py/include -I/usr/local/lib/python3.10/dist-packages/petsc4py/include -I/usr/local/dolfinx-real/include -I/usr/local/slepc/include -I/usr/local/slepc/linux-gnu-real64-32/include -I/usr/local/petsc/include -I/usr/local/petsc/linux-gnu-real64-32/include -I/usr/local/include -I/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration -I/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/wrappers -I/github/home/.cache/fenics -I/usr/include/python3.10 -c /usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/fem/DofMapRestriction.cpp -o /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/fem/DofMapRestriction.o -std=c++11 -fvisibility=hidden -std=c++20 -DDOLFINX_VERSION=0.7.0.0 -DDOLFINX_VERSION=0.7.0.0 -DHAS_ADIOS2 -DHAS_SLEPC -DHAS_PTSCOTCH INFO:root:x86_64-linux-gnu-gcc -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -fPIC -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp -I/usr/local/lib/python3.10/dist-packages/mpi4py/include -I/usr/local/lib/python3.10/dist-packages/petsc4py/include -I/usr/local/dolfinx-real/include -I/usr/local/slepc/include -I/usr/local/slepc/linux-gnu-real64-32/include -I/usr/local/petsc/include -I/usr/local/petsc/linux-gnu-real64-32/include -I/usr/local/include -I/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration -I/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/wrappers -I/github/home/.cache/fenics -I/usr/include/python3.10 -c /usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/fem/petsc.cpp -o /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/fem/petsc.o -std=c++11 -fvisibility=hidden -std=c++20 -DDOLFINX_VERSION=0.7.0.0 -DDOLFINX_VERSION=0.7.0.0 -DHAS_ADIOS2 -DHAS_SLEPC -DHAS_PTSCOTCH INFO:root:x86_64-linux-gnu-gcc -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -fPIC -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp -I/usr/local/lib/python3.10/dist-packages/mpi4py/include -I/usr/local/lib/python3.10/dist-packages/petsc4py/include -I/usr/local/dolfinx-real/include -I/usr/local/slepc/include -I/usr/local/slepc/linux-gnu-real64-32/include -I/usr/local/petsc/include -I/usr/local/petsc/linux-gnu-real64-32/include -I/usr/local/include -I/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration -I/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/wrappers -I/github/home/.cache/fenics -I/usr/include/python3.10 -c /usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/fem/sparsitybuild.cpp -o /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/fem/sparsitybuild.o -std=c++11 -fvisibility=hidden -std=c++20 -DDOLFINX_VERSION=0.7.0.0 -DDOLFINX_VERSION=0.7.0.0 -DHAS_ADIOS2 -DHAS_SLEPC -DHAS_PTSCOTCH INFO:root:x86_64-linux-gnu-gcc -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -fPIC -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp -I/usr/local/lib/python3.10/dist-packages/mpi4py/include -I/usr/local/lib/python3.10/dist-packages/petsc4py/include -I/usr/local/dolfinx-real/include -I/usr/local/slepc/include -I/usr/local/slepc/linux-gnu-real64-32/include -I/usr/local/petsc/include -I/usr/local/petsc/linux-gnu-real64-32/include -I/usr/local/include -I/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration -I/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/wrappers -I/github/home/.cache/fenics -I/usr/include/python3.10 -c /usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/fem/utils.cpp -o /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/fem/utils.o -std=c++11 -fvisibility=hidden -std=c++20 -DDOLFINX_VERSION=0.7.0.0 -DDOLFINX_VERSION=0.7.0.0 -DHAS_ADIOS2 -DHAS_SLEPC -DHAS_PTSCOTCH INFO:root:x86_64-linux-gnu-gcc -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -fPIC -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp -I/usr/local/lib/python3.10/dist-packages/mpi4py/include -I/usr/local/lib/python3.10/dist-packages/petsc4py/include -I/usr/local/dolfinx-real/include -I/usr/local/slepc/include -I/usr/local/slepc/linux-gnu-real64-32/include -I/usr/local/petsc/include -I/usr/local/petsc/linux-gnu-real64-32/include -I/usr/local/include -I/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration -I/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/wrappers -I/github/home/.cache/fenics -I/usr/include/python3.10 -c /usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/la/petsc.cpp -o /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/la/petsc.o -std=c++11 -fvisibility=hidden -std=c++20 -DDOLFINX_VERSION=0.7.0.0 -DDOLFINX_VERSION=0.7.0.0 -DHAS_ADIOS2 -DHAS_SLEPC -DHAS_PTSCOTCH INFO:root:x86_64-linux-gnu-gcc -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -fPIC -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp -I/usr/local/lib/python3.10/dist-packages/mpi4py/include -I/usr/local/lib/python3.10/dist-packages/petsc4py/include -I/usr/local/dolfinx-real/include -I/usr/local/slepc/include -I/usr/local/slepc/linux-gnu-real64-32/include -I/usr/local/petsc/include -I/usr/local/petsc/linux-gnu-real64-32/include -I/usr/local/include -I/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration -I/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/wrappers -I/github/home/.cache/fenics -I/usr/include/python3.10 -c /usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/wrappers/fem.cpp -o /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/wrappers/fem.o -std=c++11 -fvisibility=hidden -std=c++20 -DDOLFINX_VERSION=0.7.0.0 -DDOLFINX_VERSION=0.7.0.0 -DHAS_ADIOS2 -DHAS_SLEPC -DHAS_PTSCOTCH INFO:root:x86_64-linux-gnu-gcc -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -g -fwrapv -O2 -fPIC -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/pybind11/include -I/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp -I/usr/local/lib/python3.10/dist-packages/mpi4py/include -I/usr/local/lib/python3.10/dist-packages/petsc4py/include -I/usr/local/dolfinx-real/include -I/usr/local/slepc/include -I/usr/local/slepc/linux-gnu-real64-32/include -I/usr/local/petsc/include -I/usr/local/petsc/linux-gnu-real64-32/include -I/usr/local/include -I/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration -I/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/wrappers -I/github/home/.cache/fenics -I/usr/include/python3.10 -c /usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/wrappers/la.cpp -o /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/wrappers/la.o -std=c++11 -fvisibility=hidden -std=c++20 -DDOLFINX_VERSION=0.7.0.0 -DDOLFINX_VERSION=0.7.0.0 -DHAS_ADIOS2 -DHAS_SLEPC -DHAS_PTSCOTCH INFO:root:x86_64-linux-gnu-g++ -shared -Wl,-O1 -Wl,-Bsymbolic-functions -Wl,-Bsymbolic-functions -g -fwrapv -O2 /tmp/tmpsq_50hk9/github/home/.cache/fenics/.rendered.multiphenicsx_85a06d4525c3410c54a93869d159d912.o /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/fem/DofMapRestriction.o /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/fem/petsc.o /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/fem/sparsitybuild.o /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/fem/utils.o /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/la/petsc.o /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/wrappers/fem.o /tmp/tmpsq_50hk9/usr/local/lib/python3.10/dist-packages/multiphenicsx/cpp/multiphenicsx/wrappers/la.o -L/usr/local/slepc/linux-gnu-real64-32/lib -L/usr/local/petsc/linux-gnu-real64-32/lib -L/usr/local/lib -L/usr/local/dolfinx-real/lib -L/usr/lib/x86_64-linux-gnu -lslepc -lpetsc -lmpi -lmpicxx -lboost_timer -ldolfinx -o /tmp/tmpsq_50hk9/multiphenicsx_85a06d4525c3410c54a93869d159d912.cpython-310-x86_64-linux-gnu.so INFO:root:copying /tmp/tmpsq_50hk9/multiphenicsx_85a06d4525c3410c54a93869d159d912.cpython-310-x86_64-linux-gnu.so -> /github/home/.cache/fenics
nu = 0.01
def u_in_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[ # type: ignore[no-any-unimported]
petsc4py.PETSc.ScalarType]:
"""Return the flat velocity profile at the inlet."""
values = np.zeros((2, x.shape[1]))
values[0, :] = 1.0
return values
def u_wall_eval(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[ # type: ignore[no-any-unimported]
petsc4py.PETSc.ScalarType]:
"""Return the zero velocity at the wall."""
return np.zeros((2, x.shape[1]))
pre_step_length = 4.
after_step_length = 14.
pre_step_height = 3.
after_step_height = 5.
lcar = 1. / 5.
gmsh.initialize()
gmsh.model.add("mesh")
p0 = gmsh.model.geo.addPoint(0.0, after_step_height - pre_step_height, 0.0, lcar)
p1 = gmsh.model.geo.addPoint(pre_step_length, after_step_height - pre_step_height, 0.0, lcar)
p2 = gmsh.model.geo.addPoint(pre_step_length, 0.0, 0.0, lcar)
p3 = gmsh.model.geo.addPoint(pre_step_length + after_step_length, 0.0, 0.0, lcar)
p4 = gmsh.model.geo.addPoint(pre_step_length + after_step_length, after_step_height, 0.0, lcar)
p5 = gmsh.model.geo.addPoint(0.0, after_step_height, 0.0, lcar)
l0 = gmsh.model.geo.addLine(p0, p1)
l1 = gmsh.model.geo.addLine(p1, p2)
l2 = gmsh.model.geo.addLine(p2, p3)
l3 = gmsh.model.geo.addLine(p3, p4)
l4 = gmsh.model.geo.addLine(p4, p5)
l5 = gmsh.model.geo.addLine(p5, p0)
line_loop = gmsh.model.geo.addCurveLoop([l0, l1, l2, l3, l4, l5])
domain = gmsh.model.geo.addPlaneSurface([line_loop])
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [l5], 1)
gmsh.model.addPhysicalGroup(1, [l0, l1, l2, l4], 2)
gmsh.model.addPhysicalGroup(2, [domain], 0)
gmsh.model.mesh.generate(2)
Info : Meshing 1D... Info : [ 0%] Meshing curve 1 (Line) Info : [ 20%] Meshing curve 2 (Line) Info : [ 40%] Meshing curve 3 (Line) Info : [ 50%] Meshing curve 4 (Line) Info : [ 70%] Meshing curve 5 (Line) Info : [ 90%] Meshing curve 6 (Line) Info : Done meshing 1D (Wall 0.000607696s, CPU 0.000657s) Info : Meshing 2D... Info : Meshing surface 1 (Plane, Frontal-Delaunay) Info : Done meshing 2D (Wall 0.0667191s, CPU 0.065888s) Info : 2520 nodes 5044 elements
mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(
gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
gmsh.finalize()
boundaries_1 = boundaries.indices[boundaries.values == 1]
boundaries_2 = boundaries.indices[boundaries.values == 2]
viskex.dolfinx.plot_mesh(mesh)
viskex.dolfinx.plot_mesh_tags(mesh, boundaries, "boundaries")
viskex.dolfinx.plot_mesh_entities(mesh, mesh.topology.dim - 1, "boundaries_1", boundaries_1)
viskex.dolfinx.plot_mesh_entities(mesh, mesh.topology.dim - 1, "boundaries_2", boundaries_2)
V_element = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 2)
Q_element = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
def run_monolithic() -> dolfinx.fem.Function:
"""Run standard FEniCSx formulation using a mixed function space."""
# Function spaces
W_element = ufl.MixedElement(V_element, Q_element)
W = dolfinx.fem.FunctionSpace(mesh, W_element)
# Test and trial functions: monolithic
vq = ufl.TestFunction(W)
(v, q) = ufl.split(vq)
dup = ufl.TrialFunction(W)
up = dolfinx.fem.Function(W)
(u, p) = ufl.split(up)
# Variational forms
F = (nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
+ ufl.inner(ufl.grad(u) * u, v) * ufl.dx
- ufl.inner(p, ufl.div(v)) * ufl.dx
+ ufl.inner(ufl.div(u), q) * ufl.dx)
J = ufl.derivative(F, up, dup)
# Boundary conditions
u_in = dolfinx.fem.Function(W.sub(0).collapse()[0])
u_in.interpolate(u_in_eval)
u_wall = dolfinx.fem.Function(W.sub(0).collapse()[0])
u_wall.interpolate(u_wall_eval)
bdofs_V_1 = dolfinx.fem.locate_dofs_topological(
(W.sub(0), W.sub(0).collapse()[0]), mesh.topology.dim - 1, boundaries_1)
bdofs_V_2 = dolfinx.fem.locate_dofs_topological(
(W.sub(0), W.sub(0).collapse()[0]), mesh.topology.dim - 1, boundaries_2)
inlet_bc = dolfinx.fem.dirichletbc(u_in, bdofs_V_1, W.sub(0))
wall_bc = dolfinx.fem.dirichletbc(u_wall, bdofs_V_2, W.sub(0))
bc = [inlet_bc, wall_bc]
# Class for interfacing with SNES
class NavierStokesProblem(object):
"""Define a nonlinear problem, interfacing with SNES."""
def __init__( # type: ignore[no-any-unimported]
self, F: ufl.Form, J: ufl.Form, solution: dolfinx.fem.Function,
bcs: typing.List[dolfinx.fem.DirichletBC], P: typing.Optional[ufl.Form] = None
) -> None:
self._F = dolfinx.fem.form(F)
self._J = dolfinx.fem.form(J)
self._obj_vec = dolfinx.fem.petsc.create_vector(self._F)
self._solution = solution
self._bcs = bcs
self._P = P
def create_snes_solution(self) -> petsc4py.PETSc.Vec: # type: ignore[no-any-unimported]
"""
Create a petsc4py.PETSc.Vec to be passed to petsc4py.PETSc.SNES.solve.
The returned vector will be initialized with the initial guess provided in `self._solution`.
"""
x = self._solution.vector.copy()
with x.localForm() as _x, self._solution.vector.localForm() as _solution:
_x[:] = _solution
return x
def update_solution(self, x: petsc4py.PETSc.Vec) -> None: # type: ignore[no-any-unimported]
"""Update `self._solution` with data in `x`."""
x.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
with x.localForm() as _x, self._solution.vector.localForm() as _solution:
_solution[:] = _x
def obj( # type: ignore[no-any-unimported]
self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec
) -> np.float64:
"""Compute the norm of the residual."""
self.F(snes, x, self._obj_vec)
return self._obj_vec.norm() # type: ignore[no-any-return]
def F( # type: ignore[no-any-unimported]
self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec, F_vec: petsc4py.PETSc.Vec
) -> None:
"""Assemble the residual."""
self.update_solution(x)
with F_vec.localForm() as F_vec_local:
F_vec_local.set(0.0)
dolfinx.fem.petsc.assemble_vector(F_vec, self._F)
dolfinx.fem.apply_lifting(F_vec, [self._J], [self._bcs], x0=[x], scale=-1.0)
F_vec.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
dolfinx.fem.set_bc(F_vec, self._bcs, x, -1.0)
def J( # type: ignore[no-any-unimported]
self, snes: petsc4py.PETSc.SNES, x: petsc4py.PETSc.Vec, J_mat: petsc4py.PETSc.Mat,
P_mat: petsc4py.PETSc.Mat
) -> None:
"""Assemble the jacobian."""
J_mat.zeroEntries()
dolfinx.fem.petsc.assemble_matrix( # type: ignore[misc]
J_mat, self._J, self._bcs, diagonal=1.0) # type: ignore[arg-type]
J_mat.assemble()
if self._P is not None:
P_mat.zeroEntries()
dolfinx.fem.petsc.assemble_matrix( # type: ignore[misc]
P_mat, self._P, self._bcs, diagonal=1.0) # type: ignore[arg-type]
P_mat.assemble()
# Create problem
problem = NavierStokesProblem(F, J, up, bc)
F_vec = dolfinx.fem.petsc.create_vector(problem._F)
J_mat = dolfinx.fem.petsc.create_matrix(problem._J)
# Solve
snes = petsc4py.PETSc.SNES().create(mesh.comm)
snes.setTolerances(max_it=20)
snes.getKSP().setType("preonly")
snes.getKSP().getPC().setType("lu")
snes.getKSP().getPC().setFactorSolverType("mumps")
snes.setObjective(problem.obj)
snes.setFunction(problem.F, F_vec)
snes.setJacobian(problem.J, J=J_mat, P=None)
snes.setMonitor(lambda _, it, residual: print(it, residual))
up_copy = problem.create_snes_solution()
snes.solve(None, up_copy)
problem.update_solution(up_copy) # TODO can this be safely removed?
up_copy.destroy()
snes.destroy()
return up
up_m = run_monolithic()
(u_m, p_m) = (up_m.sub(0).collapse(), up_m.sub(1).collapse())
0 5.423530856245577 1 0.12945659318075842 2 0.04875205308048891 3 0.012455097181191853 4 0.010454482183770904 5 0.0008411904377640517 6 4.2733147396805684e-05 7 5.832046606630726e-08 8 2.8939001202974883e-13
viskex.dolfinx.plot_vector_field(u_m, "u")
viskex.dolfinx.plot_vector_field(u_m, "u", glyph_factor=1.0)